# this script plots the evolution of a spike in the Prinz model.

from Prinz_classes import *
from cell_test_functions import *
from pylab import *

# read in the results of Prinz's own executable
V_trace = []
t = []
for lines in open('model_test.dat').readlines():
    tokens = lines.split()
    V_trace.append(float(tokens[1]))
    t.append(float(tokens[0]))
    

# compute the results of our model
c = Ccell('test_model_cell_config')

i_time = 3500
dt = 0.05
c.exercise(timestep = dt, i_time = i_time)

# resample her results to the time array we used
V_trace.insert(0,V_trace[0])
V_trace = array(V_trace,dtype=float32) * 1000 # get her units into mV
t.insert(0,0.0)
old_t = array(t,dtype=float32) * 1000 # get her times into milliseconds

 
t = arange(0.0,old_t[-1],dt)
V_trace = interp1d(old_t,V_trace)(t)


#------------------ shift our data to get it lined up with hers ------------
spikes = spike_detect(c.V_trace())
bursts = burst_detect(spikes,intraburst_variance=5)
b_index = bursts[4][0]

her_spikes = spike_detect(V_trace)
her_bursts = burst_detect(her_spikes,intraburst_variance=5)
her_b_index = her_bursts[4][0] 

r_shift = her_b_index - b_index


def plot_range(start_time,end_time,figure_number):
    if end_time > c.total_time:
        i_time = end_time - c.total_time
        c.exercise(X_init = c.X[-1],timestep = dt, i_time = i_time)


    start_i = start_time/c.dt
    end_i = end_time/c.dt

    figure(figure_number)
    
    # construct the title of the plot
    title_string = 'Maximal conductances in mS/cm^2 \n '
    item_count = 0
    for ID in c.set_channels.keys():
        prefix = c.ID_has_prefix(ID)
        name = prefix[2:] + ": "
        conductance = c.conductances[c.channels.index(ID)]
        title_string = title_string + name + str(conductance) + '  '    
    title_string = title_string + 'leak: ' + str(c.leak_conductance)

    subplot(211)
    title(title_string,fontsize=12)
    pt = c.t[start_i:end_i]
    cv = c.V_trace()[start_i-r_shift:end_i-r_shift]
    v = V_trace[start_i:end_i]
    plot(pt,cv,pt,v)
    ylabel('Voltage (mV)')
    setp(gca(),xticks=[])
    yticks(fontsize=10)
    #subplot(312)
    #plot(c.t[start_i:end_i],c.Ca_trace()[start_i:end_i])
    #setp(gca(),xticks=[])
    subplot(212)
    # plot the difference of the two voltages (abs) on a semi-logscale
    semilogy(pt,abs(cv-v))
    # only show the xticks on the bottom plot
    yticks(fontsize=10)
    ylabel('Absolute Difference (mV)')
    xlabel('time (ms)')
    title('Mean Absolute Difference = '+str(mean(abs(cv-v))), fontsize=12 )

figure_count = 0
plot_range(2300,2800,figure_count)
